Near-linear basis set dependency handling and level shifting for RHF/UHF/DFT/UDFT by vtripath65 · Pull Request #455 · merzlab/QUICK

vtripath65 · 2026-03-09T23:26:19Z

Summary

This PR adds robust handling of near-linear basis set dependencies and level shifting to QUICK's SCF and USCF solvers, and fixes several related correctness and stability issues.

Near-linear dependency (canonical orthogonalization)

When the smallest eigenvalue of the overlap matrix falls below the OVCUT threshold (default 1e-5), QUICK now detects the near-linear dependency in fullx and sets NBSuse < nbasis. The transformation matrix x is then rectangular (nbasis × NBSuse) rather than square.

All downstream SCF operations are updated to respect this reduced dimension:

quick_overlap_module: fullx now detects near-linear dependency by scanning overlap eigenvalues against OVCUT, logs the number of linearly independent functions, condition number, and smallest eigenvalue, and builds a rectangular x(nbasis, NBSuse) via tmpU(nbasis, NBSuse) and tmpS(NBSuse, NBSuse) scratch arrays.
quick_scf_module (RHF/DFT): DIIS error vectors stored as allerror(NBSuse, NBSuse, maxdiis); DIIS transformations use hold3(nbasis, NBSuse) and hold4(NBSuse, NBSuse) rectangular intermediates; diagonalization, level shifting, and MO back-transform all operate in the NBSuse-dimensional subspace. MPI broadcasts use the correct nbasis × NBSuse and NBSuse sizes for co and E.
quick_uscf_module (UHF/UDFT): All of the above applied independently to alpha and beta spins; MPI broadcasts corrected for co/cob (nbasis × NBSuse) and E/Eb (NBSuse).

New data fields in `quick_qm_struct_type`

oeff(NBSuse, NBSuse) — reduced-space alpha Fock matrix for diagonalization and level shifting
oeffb(NBSuse, NBSuse) — same for beta spin (UHF/UDFT)
oldvec(NBSuse, NBSuse) — alpha MO coefficient matrix in the reduced space, used to apply the level shift from the previous iteration
oldvecb(NBSuse, NBSuse) — same for beta spin (UHF/UDFT)
co(nbasis, NBSuse) and cob(nbasis, NBSuse) — alpha/beta MO coefficient matrices resized from (nbasis, nbasis) to (nbasis, NBSuse)
E(NBSuse) and Eb(NBSuse) — orbital energy vectors resized from nbasis to NBSuse

A new allocate_quick_qm_struct_fullx subroutine allocates all dimension-dependent fields (called from fullx after NBSuse is determined).

Level shifting

Level shifting is applied to improve convergence of difficult cases. Three new keywords control the behavior (all can be set in the input file):

Keyword	Default	Description
`LSHIFT_CYCLE`	3	Earliest SCF cycle at which level shifting is attempted
`LSHIFT_ERR`	0.1	Minimum DIIS error required to trigger level shifting
`LSHIFT_GAP`	0.2	Desired HOMO–LUMO gap enforced by the shift

For UHF/UDFT, alpha and beta shifts are applied independently based on their respective HOMO indices (nelec and nelecb).

Overlap matrix cutoffs

OVCUT (default 1e-5): eigenvalue threshold for near-linear dependency detection.
ovmatelems (default 1e-6): elements of the overlap matrix smaller than this value are set to zero before diagonalization.

New regression tests

Four new tests covering near-linear dependency cases (all have NBSuse < nbasis):

Input	Method	Basis
`ene_benzene_b3lyp_aug-cc-pvdz`	B3LYP	aug-cc-pVDZ
`ene_benzene_ub3lyp_aug-cc-pvdz`	UB3LYP	aug-cc-pVDZ
`ene_caffeine_b3lyp_631+gd`	B3LYP	6-31+G(d)
`ene_caffeine_ub3lyp_631+gd`	UB3LYP	6-31+G(d)

Closes #436.

Remove the eigen vectors of overlap matrix corresponding low eigen function (less than 1E-05) while forming S^{-1/2}

The debug prints are removed and thresholds are returned to original values. Arrays are properly allocated to avoid multiple reallocation and memory leak.

…ifting

…lication and diagonalization

In the symmetric orthogonalization case (NBSuse == nbasis), oeff is not allocated, but the SCF DIIS procedure (electdiis) was writing to it in step 8 and the level-shift block, causing a segfault. - Step 8 transform: route the X^T.(O.X) result into o instead of oeff when NBSuse == nbasis; correct the leading dimension from NBSuse to nbasis to match o's declared shape - Level-shift block: split into canonical/symmetric branches so each operates on its own matrix (oeff vs o) - MAT_DIAG call: conditionally diagonalize oeff or o depending on the orthogonalization path - Ring-buffer, B-matrix, and C=XC' steps likewise split into canonical/symmetric branches, reusing hold2/hold instead of hold4/hold5 - Remove unused itererror staging array and hold5 scratch array - Fix sign of DGEMM 4 (alpha=-1, beta=1) to correctly accumulate e(i) = ODS - SDO without a separate subtraction loop - Add NBSuse declaration to allocate_quick_qm_struct in quick_calculated_module.f90

…w resolved

…f_module

…ole.f90

Copilot

Copilot encountered an error and was unable to review this pull request. You can try again by re-requesting a review.

agoetz · 2026-03-10T23:45:41Z

Currently, CUDA tests stop when running a test that requires f functions and the code was not compiled with f functions. This is also the case with the current master (6369ecd). I am investigating.

agoetz

This looks great. I think we need to fix minor things and Mulliken charges for two tests seem off, please see comments I made for source code lines.

src/gpu/cuda/cusolver/quick_cusolver.c

src/modules/quick_calculated_module.f90

agoetz · 2026-03-11T01:52:19Z

src/modules/quick_molspec_module.f90


      ! basis set number
      integer, pointer:: nbasis
+      integer, pointer:: NBSuse


For some reason quick_molspec%NBSuse does not really seem to be used.

agoetz · 2026-03-11T01:54:24Z

src/modules/quick_overlap_module.f90

            enddo
         enddo
-         quick_qm_struct%s(Jbas,Ibas) = SJI
+         ! do not consider small overlap matrix elements


What is the reason for not considering small overlap matrix elements?

I printed the overlap matrix elements in Gaussian and they were ignoring the small overlap matrix elements. I was experimenting with this to check if we are introducing unnecesdsary noise. I added this during my benzene aug-cc-pvdz tests which I never managed to converge with OVCUT=1.0e-6. We can remove this. Does not affect anything I think.

OK. If it does not improve convergence then we should not touch the matrix elements. I am trying to think where this could be useful or have an effect - the only thing that comes to mind is sparse matrix operations in linear scaling models. We could leave it in the code but set the default OVCUT to zero.

agoetz · 2026-03-11T02:05:27Z

src/modules/quick_scf_module.f90

    COEFF       = 0.0d0
    RHS         = 0.0d0
-    allerror    = 0.0d0
+     if(allocated(allerror)) allerror    = 0.0d0


Isn't allerror always allocated, which makes this if check superfluous?

agoetz · 2026-03-11T03:14:57Z

src/modules/quick_scf_module.f90

+                  call MAT_DGEMM ('t', 'n', NBSuse, NBSuse, NBSuse, 1.0d0, quick_qm_struct%oldvec, &
+                       NBSuse, quick_scratch%hold2, NBSuse, 0.0d0, quick_qm_struct%o, NBSuse)
+
+                  homo = quick_molspec%nelec/2


This code here is independent of dimension of NBSuse vs nbasis. Should be moved out of the if block

agoetz · 2026-03-11T03:22:03Z

src/modules/quick_scf_module.f90

+            ! Standard case (NBSuse == nbasis): diagonalize o(nbasis,nbasis) in place.
+            RECORD_TIME(timer_begin%TDiag)
+
+            if(NBSuse .ne. nbasis) then


This is where pointers can become useful. If we set

if (NBSuse .ne. nbasis) then
operator => quick_qm_struct%oeff
else
operator => quick_qm_struct%o
end if

then we don't need the if statement here, and probably also not in other places.

agoetz · 2026-03-11T03:23:03Z

src/modules/quick_scf_module.f90

-
+           ! Near-linear dependency: use hold4(NBSuse,NBSuse) as intermediate.
+           ! Standard case: use hold2(nbasis,nbasis) as intermediate.
+           if(idiis .ge. quick_method%LShift_cycle .and. errormax .gt. quick_method%LShift_err)then


if (LShift) then
...

test/saved/grad_CH4_b3lyp_def2svpd.out

test/saved/grad_SiH4_b3lyp_def2svpd.out

vtripath65 added 30 commits January 12, 2026 17:52

near-linear dependency of basis functions reduced

c36da0c

Remove the eigen vectors of overlap matrix corresponding low eigen function (less than 1E-05) while forming S^{-1/2}

simplified the handling of low overlap matrix eig values]

598cd18

rebasing after fixing cuda_diag

c710ae3

overlap matrix threshold tightened

dc5aae6

reverted the diagmkl routine back to state in master

70d5371

removed the temporary basis sets

44a436a

removed the remaining temporary basis sets

d4d2ffc

Code cleanup and proper array allocation

3d15608

The debug prints are removed and thresholds are returned to original values. Arrays are properly allocated to avoid multiple reallocation and memory leak.

MPI works

f2c6a14

CUDA code works for cases with no near-linear dependency and level sh…

922985e

…ifting

changes from master branch merged after simplifying the matrix multip…

abd18be

…lication and diagonalization

reduced the basis-set near-linear dependency default cutoff

f3467e4

use canoninical orthogonalization only in case of near-linear dependency

b94505d

some allocations are made conditional on the value of NBSuse

157bf90

opencode generated AGENTS.md

8662364

fixing compilation error due to long lines

28dd79e

missing line continuation

de673a5

The overlap matrix eigenvalue information is removed from testing

32ee544

Level shifting not applied for first two steps by default

5ed1fb1

changed some comments better readability

cf4538f

tests are added and Level shifting printing fixed

9aed47e

benzene test added to test near-linear dependency of basis sets

ecde073

Caffeine test added to test near-linear dependency of basis sets

3b3c546

Level shift keywords are added

e1c2947

the mulliken charge computation is now using MAT_DGEMM

ab944d1

the Lowdin population analysis now uses MAT_DGEMM

1021fd0

the issue with Lowdin charges in case of near-linear dependency is no…

aa8ffb4

…w resolved

fixed the test for Lowdin charges for near-linear basis set dependency

4bec1ba

Level shifting and near-linear dependency handling added to quick_usc…

3a52cfb

…f_module

vtripath65 added 6 commits March 9, 2026 15:43

added missing initializations before Lowdin charge calculation in dip…

0f8bdf6

…ole.f90

tests are added for unrestricted case

963265a

missing input test files

439ecd1

updated branch by merging master

d211b5e

removed commented out lines

3bf6db3

remove AGENTS.md (not intended for upstream)

834e1b8

Copilot AI review requested due to automatic review settings March 9, 2026 23:26

Copilot started reviewing on behalf of vtripath65 March 9, 2026 23:27 View session

Copilot AI reviewed Mar 9, 2026

View reviewed changes

agoetz self-requested a review March 10, 2026 22:23

agoetz added the enhancement New feature or request label Mar 10, 2026

agoetz requested changes Mar 11, 2026

View reviewed changes

agoetz and others added 3 commits March 10, 2026 21:29

Merge remote-tracking branch 'upstream/master' into linear_dep

52cf9cd

removed unnecessary semicolon

2eda400

Merge branch 'linear_dep' of github.com:vtripath65/QUICK into linear_dep

12d047c

Conversation

vtripath65 commented Mar 9, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary

Near-linear dependency (canonical orthogonalization)

New data fields in quick_qm_struct_type

Level shifting

Overlap matrix cutoffs

New regression tests

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Uh oh!

agoetz commented Mar 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

agoetz left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

vtripath65 commented Mar 9, 2026 •

edited

Loading

New data fields in `quick_qm_struct_type`

agoetz commented Mar 10, 2026 •

edited

Loading